import java.text.DecimalFormat;
import java.util.HashMap;
import java.util.Map;
import java.util.Random;
import java.text.*;
import java.util.*;
import java.io.*;

public class SecondOrderMM extends MM {

  protected double[][][] loge;

  private static double[][] toFirstOrder ( double[][][] emat ) {
        double[][] aux = new double[emat.length][emat.length];
        for ( int i = 0 ; i < emat.length; i++ ) for ( int j = 0 ; j < emat.length ; j++ ) aux[i][j] = 0.0;
        for ( int i = 0 ; i < emat.length; i++ ) for ( int j = 0 ; j < emat.length ; j++ ) for ( int k = 0 ; k < emat.length ; k++ ) aux[j][k] += emat[i][j][k];
        return aux;
  }

  public SecondOrderMM (Object[] esym, double[][][] emat ) {
    super(esym,toFirstOrder(emat));
    for (int i=0; i<emat.length; i++) if (esym.length != emat[i].length) throw new IllegalArgumentException("MM: esym and emat disagree");
    this.esym = esym;
    nesym = esym.length;
    loge = new double[emat.length + 1][emat.length + 1][esym.length];
    for (int b=0; b < nesym; b++) {
      double fromstart = Math.log(1.0 / esym.length );
      loge[0][0][b] = fromstart;
      for (int j=0; j < emat.length; j++) for (int k=0; k < emat.length; k++) {
	loge[j+1][0][b] = Double.NEGATIVE_INFINITY;
	loge[0][k+1][b] += emat[j][k][b];
	loge[j+1][k+1][b] = Math.log(emat[j][k][b]);
      }
      for (int k=0; k < emat.length; k++) loge[0][k+1][b] = Math.log(loge[0][k+1][b]); 
    }
  }

  public void print(PrintStream out) { 
    out.println("Symbol generation probabilities:");
    for (int b=0; b<nesym; b++) out.print(esym[b] + hdrpad);
    out.println();
    for (int i=1; i < loge.length; i++) for (int j=1; j < loge.length; j++) {
      for (int b=0; b<nesym; b++) out.print(fmtlog(loge[i][j][b]));
      out.println();
    }
  }

  public static SecondOrderMM train(Object[][] xs, Object[] esym ) {
    int nseqs  = xs.length;
    int nesym  = esym.length;
    double[][][] E = new double[nesym][nesym][nesym];
    double[][][] emat = new double[nesym][nesym][nesym];
    for (int i=0; i<nesym ; i++) for (int j=0; j<nesym ; j++) for (int k=0; k<nesym ; k++) E[i][j][k] = 0.0;
    SecondOrderMM mm = new SecondOrderMM(esym,E);
    for (int s=0; s < nseqs; s++) {
        Object[] x = xs[s];
        for (int i=0; i < x.length; i++) for (int j=1; j < x.length; j++) for (int k=2; k < x.length ; k++) E[mm.getSym(x[i])][mm.getSym(x[j])][mm.getSym(x[k])] += 1.0;
    }
    double count = 0;
    double countZero = 0;
    for (int j=0; j<nesym; j++) {
        double Eksum = 0;
        for (int k=0; k<nesym; k++) for (int b=0; b<nesym; b++) {
		count++;
		if ( E[j][k][b] == 0.0 ) countZero++;
		Eksum += E[j][k][b];
	}
        for (int k=0; k<nesym; k++) for (int b=0; b<nesym; b++) emat[j][k][b] = ( E[j][k][b] + 1.0 ) / ( Eksum + 0.0 + ( nesym * nesym ) );
    }
    if ( debugTrain ) System.err.println("Model sparsness : " + countZero + " in " + count + " ( " + ( countZero / count ) + " % )");
    return new SecondOrderMM(esym,emat);
  }

}
